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Inner magnetospheric superthermal electron transport: 
Photoelectron and plasma sheet electron sources 

/;>/■ 

G. V. Khazanov , 1 M. W. Liemohn , 2 J. U. Kozyra , 3 and T. E. Moore 4 

Abstract. Two time-dependent kinetic models of superthermal electron transport are combined to 
conduct global calculations of the nonthermal electron distribution function throughout the inner 
magnetosphere. It is shown that the energy range of validity for this combined model extends 
down to the superthermal-thermal intersection at a tew eV, allowing for the calculation of the en- 
tire distribution function and thus an accurate heating rate to the thermal plasma. Because of the 
linearity of the formulas, the source terms are separated to calculate the distributions from the 
various populations, namely photoelectrons (PEs) and plasma sheet electrons (PSEs). These dis- 
tributions are discussed in detail, examining the processes responsible for their formation in the 
various regions of the inner magnetosphere. It is shown that convection, corotation, and Cou- 
lomb collisions are the dominant processes in the formation of the PE distribution function and 
that PSEs are dominated by the interplay between the drift terms. Of note is that the PEs propa- 
gate around the nightside in a narrow channel at the edge ot the plasmasphere as Coulomb colli- 
sions reduce the fluxes inside of this and convection compresses the flux tubes inward. These dis- 
tributions are then recombined to show the development of the total superthermal electron distribu- 
tion function in the inner magnetosphere and their influence on the thermal plasma. PEs usually 
dominate the dayside heating, with integral energy fluxes to the ionosphere reaching 10 10 eV cm' 2 
s' 1 in the plasmasphere, while heating from the PSEs typically does not exceed 10 g eV cm' 2 s°. 
On the nightside, the inner plasmasphere is usually unheated by superthermal electrons. A feature 
of these combined spectra is that the distribution often has upward slopes with energy, particularly 
at the crossover from PE to PSE dominance, indicating that instabilities are possible. 


1. Introduction 

A distinctive feature of inner magnetospheric plasma is the 
presence of nonthermal electrons in the energy range of sev- 
eral eV to several keV. This electron population is formed in 
the ionosphere as a result of ionization of the atmospheric 
neutral atoms and molecules by photoionization or impact 
ionization and also convects in from the tail through the 
plasma sheet. These superthermal electrons play a very impor- 
tant role in a large number of ionospheric and plasmaspheric 
processes. 

Superthermal electrons created in the lower ionosphere 
(below about 250 km) deposit their energy before they have a 
chance to move out of this region. This is due to the high den- 
sity of neutral particles reducing the mean tree path. Inelastic 
scattering processes, combined with the source spectrum, de- 
termine the fine structure of the electron distribution function. 
The distribution is also fairly isotropic. In this region, a local 
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equilibrium approximation (which omits transport) with dis- 
crete energy loss is valid [e.g., Victor et al., 1976; kisperse, 
1976]. 

In the upper ionosphere, however, transport processes be- 
come more important. A portion of the superthermal electrons 
generated in this region can escape from the atmosphere and 
travel out along magnetic field lines. Collisions are still very 
important in determining the structure of the distribution func- 
tion, with a transition from discrete to continuous (i.c.. Cou- 
lomb collisions) energy-loss mechanisms dominating the 
energy deposition of the superthermal electrons. Transport 
models that include several pitch angle grid points, such as the 
two-stream and multi-stream approaches, have been exten- 
sively used [e.g.. Banks and Nagy, 1970; Mantas and Bowhill , 
1975; Polyakov et ai. , 1976; Strickland et ai. 1976; Prather 
et al. , 1978; Stamnes, 1980; Khazanov and Gefan , 1982; Por- 
ter et al. , 1987; Lumnierzheim et al., 1989; Link , 1992], as 
well as Monte Carlo particle tracking methods [e.g., Berger et 
al. , 1970; Cicerone and Bowhill , 1971; Solomon , 1993]. 

An important aspect of ionospheric source superthermal 
electrons is their transport through the magnetosphere, partic- 
ularly pholoelectrons (PEs) along low latitude to midlatitude 
field lines through the plasmasphere. Since Coulomb colli- 
sions of superthermal electrons with thermal electrons and 
ions are rare in this region, it was thought that PEs pass 
through the plasmasphere to the conjugate ionosphere unhin- 
dered. In this case, the plasmasphere is transparent and can be 
removed from the calculation. However, discrepancies were 
noticed between observations and this assumption [Galperin 
and Mulyarchik , 1966; Peterson et ai, 1977], and it was real- 
ized that scattering processes are sufficiently important to be 
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included in the calculation. The small-angle scattering of the 
PEs in the plasmasphere can change the motion of the elec- 
trons enough to trap some of them in the plasmasphere. This 
trapping is caused by the magnetic bottle set up by the inho- 
mogeneous magnetic field strength along the flux tube. Initial 
calculations were qualitative and focused on the thermal plasma 
heating due to this plasmaspheric scattering [ Sanatani and 
Hanson , 1970; Nagy and Banks , 1970], and later studies calcu- 
lated the differential transparency of the plasmasphere 
[Takahashi, 1973; Swartz et al. % 1975; Lejeune and Wormser , 
1976; Khazanov et ai, 1979b]. The PE distribution function 
in the plasmasphere was eventually calculated [Mantas et al . , 
1978; Lejeune , 1979; Polyakov et al. , 1979; Khazanov et ai, 
1992], but these were steady state calculations based on sepa- 
rate treatments of the ionosphere and plasmasphere. A spa- 
tially self-consistent but time-independent calculation was 
conducted by Khazanov et al. [1994], extending the results of 
Khazanov et ai [1992], Khazanov et al. [1993] presented a 
time-dependent plasmaspheric calculation, investigating the 
refilling and depletion rates of the PEs in the trapped zone. 
This model was then extended to include the ionospheres at 
each end of the field line for the first non steady state, spa- 
tially self-consistent calculation of PEs [Khazanov and 
Liemohn , 1995]. 

One process, however, that was not included in the study by 
Khazanov and Liemohn [1995] is magnetospheric convection. 
It was assumed that the flux tube was corotating with the Earth, 
so that cross field line drift effects could be neglected. Re- 
cently, simulations on a global scale have been carried out 
showing the influence of these cross field drift effects on the 
development of the high-energy PE distribution function 
[Khazanov et al. , 1996]. This study showed that this popula- 
tion requires many hours to reach a steady state level in the 
magnetosphere. 

Another source of superthermal electrons in the inner mag- 
netosphere is the plasma sheet. Because of the sunward flow of 
plasma from the dawn-to-dusk electric field in the magneto- 
sphere, plasma sheet electrons (PSEs) and ions are continually 
pushed toward the Earth, regulated by geomagnetic activity. 
Alfven and Falthammar [1963] calculated drift patterns for 
these particles, showing the classic tear-drop boundary be- 
tween closed and open trajectories (that is, trajectories that re- 
turn to their starting point and those that do not). This separa- 
trix is often called the Alfven boundary, and while the magne- 
tosphere is never constant enough for the steady-state assump- 
tion of this boundary to be valid, it is a very powerful tool for 
qualitatively describing the shape and motion of particles in 
near-Earth space. The motion of the PSEs in from the tail, 
around this boundary, and out on the dayside is well known and 
regularly observed at geosynchronous orbit [e.g., Deforest 
and Mcllwain , 1971; Eather et al ., 1976; Huliqvist et al . , 
1981; Arnoldy , 1986; McComas et al. , 1993; Birn et al . , 
1997]. Enhancements in geomagnetic activity cause this 
boundary to move in closer to the Earth, and there have been 
numerous observational and theoretical studies on the motion 
of electrons in this injection front toward the Earth near local 
midnight [e.g., Roederer , 1970; Barfield et al., 1977; Ejiri, 
1978; Ejiri et al., 1978, 1980; Kaye and Kivelson , 1979; 
Moore et ai, 1981; Moore and Arnoldy , 1982; Arnoldy and 
Moore , 1983; Kerns et al., 1994; Burke et al., 1995; Liemohn 
et ai, 1998]. A typical approach to modeling this population 
is to average the distribution function along the field line 
(with appropriate mapping to conserve the first adiabatic in- 


variant), w iich is valid if the bounce period is much shorter 
than any drift or col 1 isional timescale. This assumption is 
usually true for electrons greater than a few tens of electron 
volts. While the injected PSE distribution is usually isotropic 
(Birn et al. [1997] found that the TJT n ratio for PSEs at geo- 
synchronous is usually between 1.0 and 1.3), it has sometimes 
been noted that the pitch angle distribution of sub-keV elec- 
trons in the injection front form a possible source cone distri- 
bution, and that high-energy electrons can have a peak at 40° 
pitch angle or less [Moore and Arnoldy , 1982; Koons and 
Fennell, 1983; Arnoldy, 1986]. A strong increase in the in- 
tensity of rhe low-energy fluxes during injections has also 
been observed [Ejiri et ai, 1980]. One observation of injected 
PSEs showed the formation of a banded energy structure in the 
captured electrons [Burke et ai, 1995]. This event was simu- 
lated with a time-dependent, bounce-averaged model by 
Liemohn ei al. [1998], determining that the bands in energy 
are a natural consequence of captured electrons due to the super- 
corotation of high-energy electrons caused by magnetic gradi- 
ent-curvature drift. The bands only formed when the injected 
PSEs were contained inside the Aifvdn boundary for a pro- 
longed period. The increase in intensity in the sub-keV energy 
range was also shown to be a natural occurrence of the relaxa- 
tion of the lux tubes on the nightside. 

The theo'etical studies of this nonthermal electron popula- 
tion in the nner magnetosphere have so far not addressed the 
combinatior of the sources described above. The present study 
provides a quantitative description of superthermal electron 
distribution function formation on a global scale in the pres- 
ence of bot i PE and PSE source populations for various geo- 
magnetic conditions. In particular, the extension of the 
bounce-averaged energy range down to several electron volts 
is examinee, allowing for an accurate comparison of not only 
the relative intensities but also the energy deposition to the 
thermal plasma for the two populations. The observations of 
captured PSEs by CRRES [Burke et ai, 1995] will also be fur- 
ther examined with this combined distribution function. 

2. Model Description 

The mod;l for this study to calculate the superthermal elec- 
tron distribution function in the ionosphere-plasmasphere 
system is based on coupling two of our existing models; our 
field-alignet and bounce-averaged transport codes. These two 
models corrplement each other and for the first time offer a 
unique possibility for simulating superthermal electron mo- 
tion on a glubal scale. 

2.1. Coi figuration of the Field-Aligned 
and Bounce-Averaged Models 

The first model simulates superthermal electron transport 
through a flux tube in the geomagnetic field. It calculates the 
time-depencent superthermal electron distribution function, /, 
as a function of time, distance along the field line, energy, and 
pitch angle from the gyration-averaged kinetic equation 
[Khazanov and Liemohn , 1995; Liemohn et al. , 1997a]; 

4e dt ^ ds 2 (b ds E)dn dE\E) ] 

= Q+S 

where <t>=2Ef/m 2 is the superthermal electron flux; 
/3=(2/m) 1/2 =l .7x1 0" 8 eV 1/2 s/cm; t is time; j is the distance 
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along the field line; E is the particle energy; and \i is the co- 
sine of the pitch angle. The inhomogeneity of the geomag- 
netic field, B , is included, as well as other forces, such as elec- 
tric fields, in F. Q is the superthermal electron source term and 
S includes the collision integrals, representing interactions 
with thermal electrons and ions, scattering with neutral parti- 
cles, and wave-particle interactions. 

Initially, this model was developed for transport in the 
plasmasphere to investigate the refilling and depletion time- 
scales of the trapped zone [Khazanov et al. , 1993], with 
boundary fluxes needed at the ionospheric interfaces. The 
ionospheres were then incorporated [Khazanov and Liemohn , 
1995, 1998; Liemohn and Khazanov , 1995], including elastic 
and inelastic collisions with neutral particles, for a spatially 
self-consistent calculation along a flux tube. In these studies, 
the concept of plasmaspheric transparency was addressed, 
showing that this quantity is highly energy dependent and that 
a unified approach is needed. Recently, it was coupled with the 
time-dependent, field-aligned, hydrodynamic thermal plasma 
model of Gutter et al. [1995] to determine the coupling eifects 
between these populations [Liemohn et al ., 1997a, Liemohn 
and Khazanov , 1998], particularly the influence of a self-con- 
sistent field-aligned electric field. 

The second model of superthermal electron transport calcu- 
lates the distribution on a global scale throughout the sub- 
auroral magnetosphere. When the flight time along the mag- 
netic field line, T^, is very short compared to the collisional 
timescales, it is possible to average the particle fluxes along 
the field line over a magnetic mirror bounce period. By includ- 
ing particle drifts across field lines, superthermal electron 
fluxes for the entire inner magnetosphere can be determined by 
solving the bounce-averaged kinetic equation [Khazanov et 
al , 1996] : 

M+/ v 2) 

dt ' D/ dR ± \dtj dE \dildHo W T ft 

where includes the spatial directions R and (p perpendicular 
to the magnetic field and V D is the drift velocity in these coor- 
dinates, Hq is the cosine of the pitch angle at the magnetic 
equator, and <£> denotes averaging % over a bounce period 
along the field line. The spatial extent of this model encircles 
the globe at radial distances from L=1.75 out to L=6.5, provid- 
ing a complimentary solution of the kinetic equation to the 
first model throughout the subauroral magnetosphere. Cur- 
rently, the only collisional process is Coulomb interactions 
with the thermal plasma, and atmospheric precipitation is 
treated as a loss. 

This model has been used for PE and PSE applications. The 
formation of the high-energy PE distribution (E > 50 eV) was 
examined with this model [Khazanov et al. , 1996], showing 
that cross-field convection is an important process, particu- 
larly the transport of electrons through the nighlside into the 
dawn sector. This model has also been used to investigate the 
formation of the banded electron structure observed by Burke et 
al [1995] as PSEs were injected into the inner magnetosphere 
[Liemohn et al. , 1998]. That study focused on the capture of 
the electron cloud and the formation of the bands; this study 
will proceed a step further and examine the contribution of PEs 
to the observed distributions. 

2.2. Combined Global Model 

As was mentioned above, the simulation of the superlher- 
mal electron population is based on two major source regions: 


the ionosphere (PEs and the secondary electrons produced by 
them) and the magnetotail (PSEs). Because the PEs are gener- 
ated in the atmosphere, their plasmaspheric velocity space dis- 
tribution has a source cone at pitch angles that map to the 
ionospheres. The PE trapped population is built up by scatter- 
ing these electrons while they are in the plasmasphere. The 
timescale for the formation of the distribution function in the 
source cone is much faster than that for the trapped zone be- 
cause the former is controlled by field-aligned transport and 
the latter by diffusion [Khazanov et al. , 1993; Liemohn et al . , 
1997a], It was also shown by Khazanov et al. [1996] that 
cross-field drift effects are significant to the development of 
this trapped zone population. Because of these vastly different 
timescales, it is possible to separate the solution into two 
components: (1) a field-aligned calculation to determine the 
source function and the plasmaspheric source cone distribution 
and (2) a bounce-averaged calculation for the trapped zone. 

The PSE calculation can also be separated into these two 
components. As the plasma convects around the nightside, it 
moves closer to the Earth, experiencing a stronger magnetic 
field where a greater amount of the pitch angle domain is con- 
nected to the ionosphere. This part of velocity space is a loss 
cone, and particles will undergo severe energy loss upon enter- 
ing the upper atmosphere. This loss process is much faster 
than the convection or scattering of particles into the fly- 
through zone, and so the calculation can be divided as above: a 
bounce-averaged calculation for the development of the 
trapped zone fluxes and a field-aligned calculation for the loss 
cone and ionospheric regions, including the production of 
secondary electrons. 

This study combines these two separate calculations into a 
single consistent approach. The two models described above 
present a unique possibility in modeling the non steady state 
development of the superthermal electron distribution func- 
tion, providing complementary calculations for a feedback 
loop of loss cone-trapped zone distribution functions. The 
field-aligned model is used to calculate the PE distribution 
function in the plasmaspheric source cone, and these results 
are used as a boundary condition on the dayside in the bounce- 
averaged model. This allows tor the calculation ot the initial 
PE source terms in the ionosphere and their self-consistent 
propagation between the conjugate footpoints through the 
inner magnetosphere with the first model, and then using these 
flux intensities in the second model to simulate the cross-field 
drift effects on the PE distribution function. It also allows for 
the calculation of the PE influence to the nightside plasma- 
sphere. 

Conversely, the bounce-averaged model can be used to cal- 
culate a boundary condition for the field-aligned model. Here 
precipitation effects of superthermal electrons can be investi- 
gated, including thermal plasma heating and secondary produc- 
tion in the ionosphere. These field-aligned model results ot 
backscattering and secondary electron production can then be 
incorporated back into the bounce-averaged model to calculate 
the cross-field drift of these electrons. 

Khazanov et al. [1992] discussed an analytical approach to 
self-consistently coupling the superthermal electron solutions 
in the loss cone and in the trapped zone. In that study, a 
boundary flux injected into the plasmasphere from the topside 
ionosphere determined the distribution in '.he fly-through zone 
using a field-aligned transport equation, and a bounce-averaged 
kinetic equation was used in the trapped zone. The solutions of 
these two equations were matched at the pitch angle of the loss 
cone-trapped zone boundary. This solution was used to ealeu- 
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late the transparency of the plasmasphere for the superthermal 
electrons and from this the heating rates to the thermal 
plasma. These results, however, are dependent on the func- 
tional form of the assumed distribution at the ionospheric 
boundary. The model in the present study avoids this restric- 
tion because it contains a unified approach to the spatial re- 
gions by coupling the field-aligned and the bounce-averaged 
transport models. Also, this combined model is time depend- 
ent, so the development of the distribution function can be ex- 
amined. 

2.3. Numerical Implementation 

Each code maintains its own numerical implementation be- 
cause the only interface is the exchange of results at appropri- 
ate times. Below is a brief outline of the modeling procedure 
for each approach, which involves discrebzing the derivatives 
on a phase space grid and solving for/. A comprehensive de- 
scription has been given previously for both the field-aligned 
model [Khazanov et ai, 1994; Khazanov and Liemohn , 1995; 
Liemohn et ai , 1997a] and for the bounce-averaged model 
[Fok et al . , 1993; Jordanova et al. , 1996; Kha zanov et a l . , 
1996]. 

The field-aligned code includes the calculation of the distri- 
bution in the magnetosphere, where pitch-angle diffusion is 
the dominant process for scattering particles in to and out of 
the trapped zone. Therefore the numerical scheme of this model 
is a second-order pitch angle diffusion scheme at each point in 
time, energy, and distance, with fully implicit differences for 
the derivatives in these three variables. In order to decrease 
undesirable computational effects associated with approxima- 
tion errors of the derivatives d/ds and d/dfi, it is convenient to 
change variables from ( fj , s) to (ju 0 , s ) [Khazanov et ai, 
1979b, 1994], where /Jq~0 -( 1 ~fi 2 )B^B) xn Because ^ comes 
from the first adiabatic invariant, 0 is now a slowly varying 
function with ^ [Khazanov et al. , 1993]. The pitch angle grid 
is nonuniform, allowing more grid points near the edge of the 
loss cone where the distribution function rapidly changes, 
typically with 100 grid points per directional hemisphere. 
These grid points are added as the magnetic field decreases, 
starting with 5 to 10 in the low-altitude ionosphere and 
building up to resolve the distribution in the geomagnetic 
trap. Liemohn and Khazanov [1995] presented a calculation of 
energy conservation and numerical diffusion for this model, 
finding that the energy conservation is within the expected 
numerical error (1-10%, depending on grid spacing, when cal- 
culated to steady state). 

The bounce-averaged global model uses a combination of 
second-order high-resolution methods for the advection opera- 
tors and a second-order diffusion scheme for pitch angle scat- 
tering. The time step here must be chosen small enough to sat- 
isfy the Courant-Fredricks-Lewy condition at each phase space 
point, I v CFL \<\a{R, tp, E,/J 0 , t)At!Ax\< 1 , where Ax is the local 
step in variable x (R, (p , E, or //q) and a is the corresponding 
velocity. The time step is therefore usually less than 10 s for 
superthermal electron simulations, and At- 5 s was used 
throughout these runs. This code was recently modified for 
parallel computing efficiency, with a speed up of well over an 
order of magnitude compared to its performance on a worksta- 
tion [ McGuire and Liemohn , 1997], The numerical diffusion 
and dispersion of this technique has been discussed by Leveque 
[1992] and Fok [1993], with much less than 1% error in the 
results (when calculated to steady state). 


The field-aligned and bounce-averaged model use the same 
energy grid. For this study, it is a geometrically increasing 
energy step to provide a smooth transition from low to high 
energies. This also allows for many points at low-energies to 
capture the fine structure in this region, while simultaneously 
supplying reasonable number of grid points at high energies. 
The number of energy steps ranged from 65 to 75 for the simu- 
lations presented below. 

The linearized Fokker-Planck collision operator will be 
used for the Coulomb interactions in (1) and (2) [FI inton, 
1983, as described by Khazanov and Liemohn , 1995]. Because 
the simulations can be separated out by the source terms, that 
is, into distinct PE and PSE calculations, results from each 
population will be individually presented and discussed, fol- 
lowed by a discussion of the combined distribution from both 
sources. First, however, a critical issue must by addressed: the 
proper low energy limit of this model. 

3. Low-Energy Limit of the Combined Model 

Traditionally, bounce averaging of the kinetic equation is 
only applicable down to energies with a bounce period much 
smaller that the collisional timescales, typically above sev- 
eral tens of eV for electrons. However, once a “quasi-steady- 
state” flux ^evel is achieved in the fly-through zone (see sec- 
tion 3.2), a bounce-averaged approximation can be applied for 
much lower energies, even for electrons of a few eV energy. 
Phenomenologically, this extension is justified because of the 
very different timescales for distribution function development 
in the loss rone and trapped zone, and separate models are used 
to calculate the fluxes in these two regions of velocity space. 

3.1. Desired Low-Energy Limit 

To determine the desired low-energy limit, we will compare 
the superthermal electron energy spectra with that of a typical 
thermal electron distribution. Figure 1 shows pitch angle av- 
eraged fluxes for these two populations (thermal and nonther- 
mal electrons) at two altitudes along an L=4 field line, in units 
of(cm 2 eV ^ sr)' 1 . The superthermal spectra is from atypical 
photoelectron source, and the fluxes have been averaged over 
pitch angle for comparison with the thermal plasma spectra, 
which is assumed to be an isotropic Maxwellian, 
(f>reO =n QFe.x\\FE/E 0 )f[%7^mE^} xl2 , where the values of n 0 and 
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Figure 1. Omnidirectional fluxes for a Maxwellian thermal 
plasma and modeled photoelectrons at two altitudes along an 
L=4 flux tube. Thermal plasma parameters are given in Table 
1 . 
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Table 1. Determining the Desired Low-Energy Limit 



n 

E o 

H e 

&int 

E„ 


cm' 3 

K 

cm 

eV 

eV 

800 km 

3.1 xlO 4 

3000 

6.9xl0 7 

2.8 

2.3 

Eq. plane 

360 

5000 

1.4xl0 10 

3.9 

3.7 


Eq are given in Table 1. The intersection energies £,•„/ are 2.8 
eV at 800 km and 3.9 eV at the top of the field line. Note that 
the pitch angle averaging includes the trapped zone at the equa- 
torial plane, which is relatively empty compared to the source 
cone. By definition, there is no trapped zone at 800 km. 

Another method of determining the desired lower-energy 
limit of the superthermal electron energy range is by compar- 
ing the transport timescale z s to the collisional timescale T ee 
(the Maxwellization time). The energy where these two quanti- 
ties equate, E h> is where the distribution changes from 
"thermal” to "nontherrnar because the electrons are collision- 
dominated below this energy and transport-dominated above 
it. The timescales can be defined as 


He _ rt e 1 _ E 2 

v v\dn e /ds\ vn e O e Avn e 


(3) 


where H e is the scale height of the electron distribution, n e is 
the total electron density, and A=27te 4 lnA=2.6x 10* 12 cm 2 eV 2 . 
From (3), we get E^={An e H e ) X11 , For the same L = 4 Held line 
used in Figure 1, £/, is 2.3 eV at 800 km and 3.7 eV at the equa- 
torial plane. The relevant quantities for this calculation are 
given in Table 1. Knowing this, it is possible to analyze (1) 
and (2) to determine if the model is valid down to this energy 
limit. So, both methods reveal approximately the same low 
energy limit for the superthermal electron energy range. 


3.2. Calculation of the Limit 

The bounce-averaged approximation can only be used when 
the distribution is a slowly-varying function of field-aligned 
distance. During the initial stages of transport through a flux 
tube, this is clearly not the case, but once this initial front of 
particles has traversed the flux tube, the variation along the 
field line is greatly reduced [Khazanov and Liemohn , 1995; 
Liemohn et al. t 1997a]. It is in this ‘‘quasi-steady-state’’ limit 
where transient flows have settled that the validity of (2) will 
be explored. 

The spatial variation of the distribution function along the 
field line for a given equatorial pitch angle is proportional to 
the bounce path length for that pitch angle. Therefore, the 
most field-aligned variation in the trapped zone occurs at its 
interface with the fly-through zone ( a />). Also, (1) and (2) 
should yield the same distribution function at a because it is 
shared by both velocity space regions. Note that such a match 
can only be done if the variation of // 7 along the field line is 
small. Therefore the validity of matching these equations is 
connected to the applicability of the bounce-averaged equation 
at the edge of the loss cone and therefore throughout the 
trapped zone. 

As discussed in section 2.2, the fluxes in the fly-through 
zone reach steady state very quickly after the initial flows have 
passed because changes in the source are very slow compared 
to the interhemispherical traversal time. Therefore the time 


derivatives can be omitted for this analysis and the field- 
aligned kinetic equation becomes [Khazanov et al. y 1993], 


df An of 

a— = — 4- 
ds E dE 


£_ 

IE 2 B- fi 0 d/ifj 


B [ 

. ,.2\ % 1 


1 

Oi 
O 1 


(4) 


where the subscripts i and 0 refer to the ionosphere-plasma- 
sphere interface (where the loss cone is defined) and the equato- 
rial plane, respectively. Also, the bounce-averaged kinetic 
equation reduces to [e.g., Khazanov et al , 1992] 


d 

fyo 




M. 

<?jU 0 


+ 2b{n 0 )E?& 


- 0 


(5) 


where a and h are coefficients from the integral over a bounce 
period. Note that (5) can be derived from (4) if/ is a slowly 
varying function of s [Khazanov el al ., 1979a]. 

Note that while ji changes for a particle trajectory along a 
field line, it's equatorial value jM) is constant (except for 
changes due to collisional scattering). To combine (4) and (5) 
into one equation to solve for //>(^ob)* we can use the tact that 
the variation near /tob °t af/fyo is much bigger than the varia- 
tion of a or jjdfjiQ. So, assuming that at ^ 0 b ' ve can wr * te $lkrd 
and a{HQ)~a< then the pitch angle diffusion terms in (4) and (5) 
can be matched to yield 

( 6 ) 

dE 


^ L= A !k _\ ] _Bo_bdl 
M ds E B, a 


Now (6) can be solved for f h . For the sake of simplicity, the 
energy variation off in the low-energy superthermal range will 
be parameterized as an exponential, flE)<*txp{‘E/E s ) % and so (6) 
reduces to an ordinary differential equation resulting in 



Figure 2. Ratios of topside ionosphere values to magnetic 
equator values of the electron distribution function at the loss 
cone-trapped zone boundary at various L shells. Flux ratios are 
from the solution of (12), with (top) E s = 10 eV (bottom) and 20 
eV. The relevant parameters given in Table 2. A heavy dashed 
line at unity is shown for reference. 
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Table 2. Loss Cone Boundary /Ratios at 4 eV 


L 

B,IB o 

«0b 

deg 

^0b 

cm* 3 

fbitfbQ fbiV\ bO 

£ 0 =10eV £ 0 =2 OeV 

2 

8.53 

20.0 

0.940 

2.3x1 0 3 

1.7 

1 .3 

3 

32.1 

10.2 

0.984 

l.OxlO 3 

2.4 

1 .6 

4 

79.8 

6.4 

0.994 

360 

2.4 

1.6 

5 

160. 

4.5 

0.997 

180 

2.3 

1.5 

6 

281. 

3.4 

0.998 

15 

1.1 

1.1 


fb,i 

fb,o 



B,j f V(s') ) 

\ B{*')EE S 

o 

oeT 



(7) 


Plots of (7) are shown in Figure 2 for £y=l0 and 20 eV at 
several L values. These choices for E s are typical mean values 
of superthermal electrons in the low-energy range. The distri- 
bution function ratio is the value of / at the loss cone boundary 
at 800 km to that at the equatorial plane. The thermal plasma 
densities used are given in Table 2 along with the magnetic 
field ratio and the equatorial pitch angle of the loss cone- 
trapped zone boundary cr 0 b ( ar >d its cosine, ^i 0b ). Figure 2 
shows that the//, ratios are quite reasonable down to the desired 
low-energy limit (the 4 eV values listed in Table 2). They are 
well below the corresponding magnetic field ratios, indicating 
that the decrease along the field line due to collisions is rather 
insignificant compared with the changes due to the inhomoge- 
neous magnetic field. Note that the L = 3 and L = 4 results are on 
top of each other, and are slightly higher than the L - 5 ratios. 
Numerical results from the field-aligned transport model con- 
firm that ob) converges very quickly, the greatest field- 
aligned variation in the trapped zone occurs at and the 
code yields very similar tlux ratio results for times after the 
initial front has passed through the plasmasphere. It is clear 
from this analysis that a bounce-averaged approach is valid 
down to very low superthermal electron energies, as long as 
impulsive field-aligned Bows are omitted. The L shell depend- 
ence of these ratio curves is mainly due to the thermal density. 

This extension allows for the calculation of the total energy 
spectrum of the nonthermal electron population, as well as an 
accurate determination of the energy deposition rate from the 
PE and PSE populations to the thermal plasma. These topics 
will now be addressed. 


4. Photoelectron Distribution Function 
Formation 

As discussed above, the linearity of (!) and (2) affirms that 
the distribution functions from the various source terms can be 
calculated independently, and then summed at the end to show 
the combined result. The PE source will be discussed first, fol- 
lowed by the PSE source term. 

4.1. Photoelectron Source Spectrum 

An important feature of the PE source spectrum is the pres- 
ence of spikes due to specific emission lines in the solar EUV 
spectrum. A prominent line is the He II 30.4 nm line. With an 
energy of 40 eV, this strong line produces a series of peaks in 
the 20-30 eV electron energy range due to the various ioniza- 
tion states of the atmospheric constituents. 


Another mportant feature of the photoclectron source spec- 
trum is the presence of the Auger electron peaks. This is due to 
the double ionization of atmospheric particles by high-energy 
photons. These photons initially produce a free electron from 
an inner electron shell of the neutral. Then an electron from 
the outermost shell transitions down to fill the hole, releasing 
a photon equal to this transition energy. This photon, how r - 
ever, does not escape but instead produces another free electron 
from the outermost shell. This second electron is known as an 
Auger electron, and has an energy of 300-500 eV, depending 
on the particle it originated from. There have been only a few 
modeling el forts to include this population [Avakyan et at., 
1977; Winn.ngham et al., 1989], showing that the inclusion 
of Auger electrons is necessary to obtain good agreement with 
observations while maintaining reasonable soft X ray solar 
fluxes. For the present study, which focuses on the compari- 
son of PE fluxes with PSE fluxes in the inner magnetosphere, 
Auger electrons are a critical component of the results. 

It is useful to compare these model results with data to de- 
termine the accuracy of the PE source spectrum to be used in 
this study, "igure 3 shows such a comparison. Here omnidi- 
rectional electron measurements from the Atmospheric Ex- 
plorer E satellite are plotted along with a model run for the 
same geophysical conditions. The data were taken on the 
morningside near the equator on day 355 of 1975 [Doering et 
al., 1976], a:, solar zenith angles of 50“ and 37° for the two al- 
titudes of 182 km and 365 km, respectively. In Figure 3a, the 
spectra agree closely for most of the energy range. Figure 3 b 
also shows good agreement, with the model predicting more 
definition in the 20-30 eV range (from photoionization from 
the He II 30 4 nm solar line) and slightly lower fluxes above 
30 eV. Furher comparisons of results from this model with 
observation? are shown by Khazanov and Liemohn [1998], 

As discu ised above, the field-aligned model provides a 
source term for the bounce-averaged model at the loss cone 



0 20 40 60 

Energy (eV) 


Figure 3, Comparison of model results (solid lines) with 
AE-E data (dashed lines) at (a) 182 km and (b) 365 km on day 
355 of 1975. The satellite data are reproduced from Doering et 
al. [1976]. 
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Figure 4. Typical photoelectron spectrum from the field- 
aligned code used as an input condition for the bounce-aver- 
aged model at the loss cone-trapped zone boundary. 


boundary on the dayside for the PEs. For the cases to be exam- 
ined here, a moderately high solar activity level will be used 
along with various geomagnetic activity levels. The PE source 
is mostly a function of solar activity level, and largely inde- 
pendent of geomagnetic activity. Therefore the PE source term 
being supplied to the bounce-averaged model is essentially the 
same for the various simulations. A typical source cone 
boundary spectrum is shown in Figure 4. The low-energy 
range below 5 eV is sloped upward due to Coulomb collisional 
damping, the production peaks visible in Figure 3 are essen- 
tially washed out into a single increase in flux near 25 eV, and 
the Auger peaks at 360 and 500 eV are quite distinguishable. 
The production spectrum drops off rapidly above 500 eV be- 
cause of the absence of the Auger electrons. Note that there is 
actually fine structure present in the low- and high-energy pro- 
duction peaks resulting from the individual ionization states of 
the atmospheric constituents, and the resolution of this fine 
structure is highly dependent on the energy step of the model. 
The choice of a geometrically increasing energy step for these 
calculations smoothes out the peaks, especially in the high- 
energy range. Choosing a AE of much less than 1 eV is neces- 
sary to accurately resolve each peak [cf. J asperse and Smith , 
1978], and the stability of these peaks have been debated as a 
source of plasma waves [Butvin et al. , 1975; Kudryashev et al., 
1979; Ivanov et al. , 1980, 1982], Because we are focusing on 
the global formation and evolution of the distribution function 
rather than the possible instabilities, a reasonably large en- 
ergy step was chosen. 

4.2. Photoelectrons When Kp~\ 

A comprehensive calculation of the superthermal electron 
distribution function throughout the inner magnetosphere 
from the thermal energy boundary up to tens of keV energies 
has not been conducted, to our knowledge. The combined 
model presented above now makes this type of time-depend- 
ent, spatially three-dimensional simulation possible for all 
superthermal electron sources. Because the PE source term is 
largely independent of geomagnetic activity, a simulation 
with constant Kp=l will be discussed in detail, and then other 
simulations will be compared to these results. 

The development of the source cone region of velocity 
space on the dayside has been shown in previous studies [e.g., 
Liemohn and Khazanov , 1995]. Here we will focus on the for- 
mation of the trapped zone flux population. The time devel- 


opment of the 90° pitch angle (equatorially mirroring) energy 
spectra is shown in Figure 5 for three L values and six mag- 
netic local times. Note the different flux scales for each sub- 
plot. Here and throughout the remaining figures, the differen- 
tial number flux is given in units of (eV cm 2 s sr)" 1 , and the dis- 
tributions are shown in terms of equatorial pitch angle (these 
equatorial plane distributions from the bounce-averaged model 
can easily be mapped along the field line). Because the source 
is at small pitch angles, the development of the distribution at 
90° should be a maximum time scale for PEs. Note the fast 
buildup of the fluxes on the dayside. A steady state distribu- 
tion has been reached for most of the PE energy range after 3 
hours. This case has relatively high thermal plasma densities, 
accelerating the scattering of PEs into the trapped zone. On 
the nightside, the PE fluxes degrade as they move around from 
dusk to dawn. The peak near 100 eV is due to the source drop- 
ping off at higher energies and Coulomb collisions more effi- 
ciently absorbing the energy of the slower PEs. The L - 5 
fluxes are much higher than the L= 3 fluxes at MLT=0300. This 
is not only because of the decreased thermal plasma density at 
the larger L values but also because of magnetospheric convec- 
tion compressing the flux tubes from larger L shells. Thus it 
appears the flux at this L shell remains constant (and even 
grows) around the nightside. At MLT=0600, this narrow L 
shell band of PE fluxes contributes to the dayside fluxes in the 
high-energy trapped zone, creating a bulge in the distribution 



Figure 5. Time development of the 90° equatorial pitch an- 
gle energy spectra at various spatial locations in the inner 
magnetosphere. Simulation times are in hours from photoe- 
lectron source turn-on in the dayside ionosphere, with con- 
stant Kp= 1 geomagnetic activity. Note that each subplot has 
its own dynamic range. Here and elsewhere, the axis label 
"flux’' refers to differential number flux in (eV cm 2 s sr)’ 1 . 
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Figure 6. Energy dependence of the steady stale equatorial 
pitch angle distributions at various spatial locations for a pho- 
toelectron source with constant Kp- 1. Note that each subplot 
has its own dynamic range. 


function. This bulge is seen even at MLT=1200. By 24 hours 
of simulation time, the fluxes have essentially reached a 
steady-state level. 

The pitch angle distributions at the spatial locations of 
Figure 5 are shown in Figure 6 for the time-converged results. 
On the dayside, it is clear that the trapped zone for low-energy 
electrons is filled in more than the higher energies. On the 
nightside, the low-energy fluxes are nearly isotropic in the 
trapped zone. This is not the case at higher energies. In the 
220 eV and 490 eV nightside distributions, there is a peak near 
30 e -40\ This is largely due to Coulomb interactions. At 
small pitch angles, scattering depletes the loss cone and 
nearby trapped zone. At large pitch angles, scattering from 
the source cone on the dayside is not able to fill the loss cone 
any more than this, because of the competing processes of 
Coulomb energy decay and spatial drift. Please note that this 
is different from the source cone distributions observed by 
Moore and A mold y [1982], which are due to secondary elec- 
trons from precipitating plasma sheet particles, a source not 
included here. Also, the bulge from the particles conveciing 
through the nightside is clearly seen at MLT=0600, contribut- 
ing only to the high-energy trapped zone fluxes. 

4.3. Analysis of Photoelectron Transport 

Figures 5 and 6 presented the development of the PE distri- 
bution function for low geomagnetic activity with all drift 
processes included. It is useful to examine the change in these 
results under different conditions and with various processes 


omitted. This will provide insight into why the distribution 
develops as a does. 

First, the influence of the various drift terms will be exam- 
ined. Using the Volland-Stern convection model [ Volland , 
1973; Stern , 1975] and a dipole magnetic field, the bounce-av- 
eraged velocity terms in (2) from adiabatic drift can be written 
in the form [e.g., Jordanova , 1995] 

ldR\ = _AR\os(p_ ( 8a) 

\ dt 1 M e 

/d(p\ C | Aft 3 siny 3 ER ] l{v o) (8b) 

\dtl M e M e qM E [ 6 / i ( mo )_ 

ldE\_ 3 E\ /(M l (dR\ (8c) 

\ dt I R 6h(n 0 ) \dt I 

(dfip\ (l-Mo) l[fto) ldR\ (8d) 

\dt/ 4 Rfi 0 h(n 0 )\dt/ 

where R is geocentric distance in the equatorial plane, <p is 
MLT in radians (<p= 0 at MLT=0000), A is a function of geo- 
magnetic activity (Kp in this case, taken from Maynard and 
Chen [1975]), C is a constant characterizing the corotation 
electric field, M£=8.02xl0 15 Tm^ is the magnetic dipole of the 
Earth, q is the charge of the particle including sign, and I(jUq) 
and hip o) ate slowly varying functions of equatorial pitch 
angle result ng from the bounce-averaging process [Ejiri, 
1978], The various processes can be identified in the right- 
hand sides o' these equations. Corotation yields the first term 
in (8b); magnetic gradient-curvature drift produces the third 
term in (8b), and the rest of the terms in (8) are due to magne- 
tospheric convection. By systematically removing these 


Sim ^ 1 Sim ^2 Sim ^3 
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Figure 7. Comparison of the time development of the 
nightside pit :h angle distributions at various energies for 
three different simulations at MLT=0300, L= 5. All three have 
constant Ay?=U, and Sim 1 includes all drift terms, Sim 2 is 
without magnetospheric convection drift, and Sim 3 is without 
magnetic gradient-curvature drift. 
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terms from the calculation, the effects of each drift process can 
be investigated. 

Figure 7 shows a comparison of the pitch angle distribu- 
tions for three simulations at a nightside spatial location in 
the PE band (MLT=0300, L- 5). All three simulations have a 
constant Kp= 1 with various combinations of the drift terms 
(listed in Table 3). Sim 1 and 3 look very similar because the 
energies are too low to be affected by magnetic drifts. Sim 2 is 
lower than the other two because without magnetospheric con- 
vection there is no relaxation on the nightside to enhance the 
flux at this L shell. These results indicate that it is only the 
high energy fluxes that benefit the most from this flux tube 
compression; the low energy region has a large Coulomb loss 
term that offsets some of the increase. Note that at higher en- 
ergies, the peak at 30°-40° pitch angles is evident, particularly 
in the early time results. Its existence in the Sim 2 results 
suggests that Coulomb collisions are mainly responsible for 
this peak. Another alternative is the field-aligned acceleration 
due to the relaxation of the stretched flux tubes (lirst-order 
Fermi acceleration). This is perhaps partly responsible, but 
Coulomb collisions play a much larger role in forming this 
peak. The 5 eV fluxes are nearly isotropic, even here on the 
nightside. Note that the dashed line is at 9 h of simulation 
time, when the PEs should first be reaching this spatial loca- 
tion. However, most of the energy range needs at least another 
3 hours to approach a steady-state level (the high-energy range 
needs much more than this, in fact). This can be explained by 
the fact that it takes some amount of time to scatter particles 
into the trapped zone in the afternoon sector (the source region 
for PEs on the nightside), so 9 hours is premature to expect a 
complete distribution at this spatial location. Also, the total 
convection rate is less than the purely corotational convection 
rate because magnetospheric convection acts to slow the PEs 
before local midnight (MLT=1 800-2400) for a longer interval 
than it acts to speed them up (MLT=2400-0300). 

It is useful to show a comparison of the relative strengths of 
the azimuthal (MLT) drift terms from (8b). These drifts are 
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Figure 8, Radial dependence of various azimuthal (MLT) 
drift terms in hours of magnetic local time per simulation hour. 
Terms shown are for corotation (solid line), magnetospheric 
convection at MLT=0600 (where the azimuthal component is 
greatest) for Kp= 1 (dotted line) and Kp = 6 (dashed line), and the 
gradient-curvature effects from the magnetic field for electron 
energies of 200 eV (dash-dot line) and 20 keV (dash-dol-dot- 
dot line) for 80° pitch angle and also 20 keV with 20° pitch 
angle (long dash line). 


Table 3. PE Comparison Simulations 


Simulation 

Kp 

Processes 

1 

1 

all included 

2 

1 

no magnetospheric convection 

3 

1 

no gradient-curvature drift 

4 

1 

no drifts at all 

5 

6 

all included 


shown in Figure 8, in units of hours of magnetic local time per 
simulation hour. Corotation is at unity and constant with L. 
The magnetic gradient-curvature drifts are weakly pitch angle 
and L dependent and strongly energy dependent. Two energy 
values are shown for a pitch angle of 80°, with the higher one 
repeated for 20°. From this it is clear that magnetic drifts are 
unimportant at PE energies, but it will be the dominant dritt 
term at higher energies. Note that due to the negative charge 
of electrons, this drift is in the same direction as corotation. 
Drift due to magnetospheric convection is dependent on Kp. 
and a high and low value are shown. This drift is also depend- 
ent on MLT, and so MLT-0600 was chosen to show the maxi- 
mum positive azimuthal drift due to convection. This <p drift 
will vanish at MLT=0000 and 1200, and will counteract coro- 
tation at MLT=1 800. This shows that corotation always 
dominates convection in the inner part of the spatial domain, 
but during high activity convection will dominate in the outer 
part of the spatial domain. This will most likely be true even 
for a more realistic convection electric field description. 

Not shown in Figure 7 are the other two simulations listed 
in Table 3. This is because the fluxes are zero for these two 
cases at this spatial location. No PEs should he expected on 
the nightside when there is no spatial drift, as is the case for 
Sim 4. When Kp - 6 (Sim 5), corotation only dominates very 
close to the Earth, and so this spatial location is outside of the 
Alfven boundary for all PE energies and no PE flux is expected. 
Also, the corotating region where PEs are expected has a large 
thermal plasma density, and so the PE fluxes in the nightside 
plasmasphere are very small. 

The effects of these two simulations will make a difference 
on the dayside, however, and this is shown in Figure 9. Here 
steady state pitch angle distributions are presented for L - 5 at 
several energies and magnetic local times. Notice that Sim 5 
is drastically different from the other results. This is because 
the enhanced magnetospheric convection pushes the plasma 
out the dayside boundary before it can build up. Only in the af- 
ternoon sector at low energies is there enough time for Cou- 
lomb scattering to fill in the trapped zone to a level close to 
the low geomagnetic activity simulations. Sim 4 and Sim 5 do 
not show the high-energy bulge on the dawnside because there 
is no plasma entering from the nightside for these cases. Sim 
2 has a much smaller trapped zone bulge intensity than Sim 1 , 
and it does not persist as far into the dayside either. Sim 2 
does, however, have a larger flux at high energies in the after- 
noon sector because there is no magnetospheric convection to 
strip the plasma away. This is also true for the comparison of 
Sim 1 and Sim 4 at low energies at dawn. Because MLT-0600 
is at the beginning of the PE source region, all azimuthal drifts 
act to deplete the low-energy trapped zone at this MLT. Thus 
Sim 4 has slightly higher fluxes at low energies here. 

Several conclusions can be made from these comparisons. 
One is that magnetospheric convection plays a large role in 
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Plate 1. Flux tube-integrated energy deposition rates to the thermal electrons for photoelectrons, plasma 
sheet electrons, and the combined distribution for three plasma sheet densities. The first row is constant Kp= 1 
at steady state, the second row is constant Kp - 6 at steady state, and the last three rows are at times during late 
January 1991. 


the development of the PE distribution. It creates a high-en- 
ergy bulge in a narrow L band on the nightside that propagates 
into the dawn sector, and it acts to deplete the trapped zone 
fluxes on the dayside by sweeping the plasma away before it 
can build up. Magnetic gradient curvature drift is not impor- 
tant for PEs, but Coulomb collisions are very important. This 
interaction fills in the trapped zone, allowing for a PE popula- 
tion on the nightside, and it removes PEs from the inner plas- 
masphere on the nightside. High geomagnetic activity com- 
pletely changes the distribution of PEs throughout the inner 
magnetosphere, preventing them from moving to the night- 
side and removing them from the dayside as well. 


5. Injection of Plasma Sheet Electrons 

For PSEs, ihe various drift terms play an even more impor- 
tant role in tie distribution function formation in the inner 
magnetosphere because of the characteristics of the source 
population. Far this reason, a real Kp history will be used for 
the baseline case, and then other simulations will be compared 
with this one :o examine the effects of the various processes. 

5*1. Plasma Sheet Boundary Conditions 

Liemohn ei ai [1998] discussed the necessity for careful de- 
termination of where the Alfv6n boundary lies in relation to 
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Figure 9. Comparison of the steady-state photoelectron 
pitch angle distributions at various energies and dayside mag- 
netic local times at L= 5 for four different simulations. Sim 1 
and Sim 2 are the same as in Figure 9, while Sim 4 has Kp= 1 
with no drift terms included at all (only collisional processes), 
and Sim 5 has constant Kp- 6 with all processes included. 


the outer boundary of the simulation domain. The plasma trav- 
els along closed drift paths inside of the Aifven boundary and 
along open drift paths outside of it, and this separalrix is 
highly dependent on geomagnetic activity, as well as energy, 
pitch angle, and location. The location in (p ot where the 
Aifven boundary crosses the outer boundary ot the simulation 
at 6.6 R E , (p B (E, fi a , Kp), must be accurately determined, with 
injected plasma sheet particles entering along the open trajec- 
tories (dawnward of (p B ) and corotating plasma distributions 
from the dayside filling closed trajectories (duskward of (p B ). 
As in the work of Liemohn et al. [1998], (p B is taken from a 
look-up table so that the appropriate boundary flux is taken for 
each energy and pitch angle at each azimuthal angle for each 
time step. 

The models can assume any initial or boundary condition. 
For the PSEs, the injection flux is assumed to be an isotropic 
kappa distribution, which resembles a Maxwellian below the 
characteristic energy Eq but transitions to a power law ot fc*E~ K 
at higher energies. While the injected electrons sometimes 
have a characteristic energy up to a few keV [DeForest and 
Mcllwain , 1971; Eather et al ., 1976; Evans and Moore , 1979; 
Birn et al ., 1997; Borovsky et al ., 1997], Christon et al 
[1989] found that best fit values for electrons from ISEE 1 
measurements during quiet times are K= 6 and Eq- 200 eV, and 
these values will be used in the present study. The density ot 
the injected electrons is taken to be 0.1 cm' 3 [cf. Huang and 
Frank , 1986; Baumjohann et al. , 1989; Christon et al , 1989; 
Birn et al. , 1997]. At the beginning of each simulation 
(UT=0000 of 1/25/91), it is assumed that there are no PSEs in- 
side of the simulation domain. 

The geomagnetic activity level is a critical element in elec- 
tron transport in the inner magnetosphere. Thus the choice of 
Kp history is very important. Several simulations will hold 
the activity level fixed at a low (Kp- 1) or high (Kp- 6) value. 


but these are simply illustrative cases. For a realistic simula- 
tion, we will use the Kp history of late January 1991. During 
this time period, the low-energy plasma analyzer (LEPA) in- 
strument on-board the CRRES satellite detected several bands 
of magnetically trapped 0.1-30 keV electrons inside the plas- 
mapause. Liemohn et al. [1998] examined the formation of 
this banded structure in the inner magnetosphere, concluding 
that these bands are a natural occurrence when plasma sheet 
electrons are trapped along closed drift paths in the inner mag- 
netosphere by a sudden decrease of geomagnetic activity. Fur- 
ther simulation results during this time period will be pre- 
sented, and the possibility of photoeleclron additions to the 
low-energy band of the trapped population causing its mean 
energy to decrease diabatically (as discussed by Burke et al. 
[1995]) will be investigated. The Kp history for late January 
1991 is shown in Figure 10. 

Notice in this plot that there are a series of spikes in Kp on 
top of a relatively low background level. Of interest are the 
three most prominent spikes at the ends of days 25, 26, and 
27. Burke et al. [1995] stated that the lowest-energy band was 
not observed by CRRES until after the last ot these spikes. 
Liemohn et al. [1998], however, found that the cloud exists 
even after the first Kp spike and that the last spike simply adds 
to the intensity of the cloud. This difference could be from a 
number of possibilities: the satellite missed the low-energy 
cloud until the 27th (and thus only measured the high-energy 
bands elsewhere in space); scattering processes not included in 
the model rapidly depleted the captured fluxes (such as wave- 
particle interactions); or the choice of PSE boundary flux is 
not consistent with the actual injections. Birn et al. [1997] 
showed that both the density and temperature of the electrons 
during an injection typically rise and then fall ott. However, 
they fitted their distributions with a Maxwellian rather than a 
kappa distribution, while the changes in temperature are most 
likely due to enhancements in the tail ot the distribution (as 
seen in their Figure 9). A Maxwellian has a higher characteris- 
tic energy than a kappa distribution would need, and thus our 
choice of boundary condition is not unreasonable. Still, an 
increase in either K or Eq would decrease the number of elec- 
trons in the low energy range, and this could explain the ab- 
sence of the low-energy band in the CRRES observations. 

5.2. Plasma Sheet Electrons During the CRRES 
Observations 

Using the Kp history in Figure 10, the PSE injection into 
the inner magnetosphere can be simulated for the CRRES ob- 
servations of the banded electron structure event. 



Day of January 1991 

Figure 10. history during the CRRES observations of 
the banded electron structures. 
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Figure 11. Radial dependence of 90° plasma sheet electrons 
at various energies and simulation times for late January 1991, 
shown at local midnight (solid line) and noon (dotted line). 


This L shell dependence is shown in Figure 1 1, illustrating 
the radial variation of the 90° (equatorially mirroring) flux val- 
ues at several energies and simulation limes. The times shown 
are during the first and third pulses and then at the end of the 
simulation, chosen so the peak of the cloud is at local mid- 
night. It is clear that essentially no low-energy PSEs pene- 
trate inside of L=3, but a portion of the high-energy tail of the 
kappa distribution makes it in this far after the third injection. 
The 650 eV row shows the fluxes near the peak of the lowest 
energy band. The degradation of the low-energy fluxes (below 
650 eV) due to Coulomb collisions is clearly seen in this fig- 
ure, although the degradation timescale is slower than for PEs 
because of the higher mean energies of the PSEs. The high-en- 
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Figure 12. Local midnight plasma sheet electron pitch an- 
gle distributions at various energies and simulations times for 
late January 1991, shown at L=4 (solid line), 5 (dotted line), 
and 6 (dashed line). 
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Figure 13. Comparison of plasma sheet electron energy 
spectra at L—4, MLT=0000 tor three different simulations at 
various times. All three have the real Kp history, but Sim 1 
has all processes included, Sim 2 has no Coulomb collisions, 
and Sim 3 1 as no Coulomb collisions or gradient-curvature 
drift (dashed line). 


ergy fluxes :n the third column look very different from the 
others because of their faster rotation period (as discussed 
above in Figure 8), so this is not a cut through the peak of the 
distribution at this energy. 

The pitch angle distribution of the injected PSEs is shown 
in Figure 12. Note that the injected distribution was isotropic 
in the trapped zone, and these results indicate that it remains 
quite isotrop c throughout the simulation in both space and 
time. The o ily variation from isotropy is at high energies af- 
ter several di.ys of rotation. This is due to the pitch angle de- 
pendence of the gradient-curvature drift. The pitch angle dis- 
tribution near the peak of the low-energy band shows signs of 
degradation after several days of rotation, and the fluxes near 
the loss cone are about an order of magnitude lower than at 90°. 

5.3. Analysis of Plasma Sheet Electron Transport 

Some analysis has already been discussed, but here a com- 
parison of the results from the previous section with other 
simulations s presented. Figure 13 shows a comparison at 
L=4 and ML”=00 of three simulations at several times. The 
simulations a'e described in Table 4. All three have the same 
Kp history, but they have different combinations of included 
processes. Here it is seen that Coulomb collisions form the 
low-energy t£ il of the captured cloud. Without Coulomb colli- 
sions (Sim 2) there are no electrons below 200 eV at this spa- 
tial location due to the Fermi acceleration experienced during 
the Earthwarc convection. This low-energy tail forms at the 
expense of tht peak intensity level of the lowest energy band. 
This is partic llarly evident in the last subplot, as the low-en- 
ergy peak is i ot degrading in time or space. Coulomb scatter- 
ing of electro] is at energies of a few hundred eV has a timescale 


Table 4. PSE Comparison Simulations 

Simulation Description 

1 real Kp y all processes 

2 real Kp , no Coulomb collisions 

3 real Kp, no CC or G-C drifts 

4 real up to t=48 then Kp= 1 , all processes 

5 Kp - 6 always, all processes 






KH AZANOV ET AL.: SUPERTHERMAL ELECTRON TRANSPORT 


23,497 



Figure 14. Comparison of plasma sheet electron energy 
spectra at L= 4, MLT=0000 for three different Kp history simu- 
lations at various times. All three include all processes, but 
Sim 1 used the real Kp history, Sim 4 used the real Kp history 
through 1/26 and then switched to constant Kp- 1, and Sim 5 
used constant Kp= 6. 


of days [Liemohn et al. y 1997a, b], and the slow diffusion ol 
particles into the loss cone is responsible for the slight degra- 
dation of the peak seen in Figure 11, With the removal ol 
magnetic gradient-curvature drift as well, there is no banded 
structure formation at all, and Sim 3 simply shows an acceler- 
ated kappa distribution. The cloud is still captured in this 
simulation, but there is no degradation or distortion of the 
fluxes. 

Figure 14 shows a comparison ot Sim 1 with two more 
cases. Here, the three cases have all processes included but 
have different Kp histories. Sim 4 is the case without the third 
Kp spike, so all of the captured electrons are from the first two 
spikes. It is clear that this final Kp pulse only adds to the ex- 
isting cloud and is not the sole source of the captured elec- 
trons, at least for the chosen boundary distribution. The dif- 
ferences between Sim 1 and Sim 4 in the last subplot is mainly 
due to differences in the thermal plasma due to the differing Kp 
histories. Sim 5, however, is completely different. The con- 
stantly high activity level pushes in the kappa distribution, 
and then constantly replenishes it. This is different from Sim 
2 in that these particles are freshly injected instead of captured 
and undamped. There is no capture in Sim 5 because the activ- 
ity is constant; a relaxation is necessary to capture a cloud of 
PSEs. 

6. Combined Electron Distribution Function 

The combination of the two sources described above gives 
new insight into the development of the total superthermal 
electron distribution in the inner magnetosphere. Also, a new 
feature of this model is the ability to calculate the low energy 
range of the superthermal electron distribution with a bounce- 
averaged technique, allowing for an accurate calculation of the 
energy deposition to the thermal plasma on a global scale. 
These two features will be discussed in this section. 

6.1. Heating the Thermal Plasma 

When energetic plasma is streaming through a background 
thermal plasma, the energetic particles will lose energy to the 
thermal plasma through Coulomb collisions. This results in 


an energy deposition from the energetic particles to the low- 
energy populations. The instantaneous volume heating rate Q 
is derived by integrating the Fokker-Planck collision term of 
Coulomb interactions over velocity space, and is a function of 
the superthermal electron omnidirectional flux and the thermal 
plasma density. Of particular interest is the heat flux into the 
ionosphere P& obtained by integrating Q along the Field line 
from the ionospheric footpoint to the equatorial plane 
[Liemohn and Khazanov, 1995]. This represents the energy 
flux into the ionospheric thermal population from the ener- 
getic population. Studies have shown that variations in P e 
have great effects on the thermal, compositional, and optical 
structures in the ionosphere [Chandler et a/., 1988; Kozyra et 
al.y 1990; Comfort et ai , 1995]. 

Integrated thermal electron heating rates are shown in Plate 
1 (it is a straightforward calculation to produce results tor the 
other thermal plasma species as well). The five columns are 
defined as follows: heating from PEs; heating from PSEs; 

summed PE and PSE heating (sum of the first two columns); 
summed PE and PSE heating with 5 times the PSE density at the 
boundary as in the previous two columns (so /ipsE = 0>5 cm’ 3 in 
the boundary kappa distribution); and summed PE and PSE 
heating with «pgp=3.0 cm' 3 at the boundary (30 times larger). 
This increase in density is simply a multiplier because of the 
linearity of (1) and (2), assuming that the injected distribution 
is held constant throughout the simulation. These last two 
columns are shown to represent injection from a dense plasma 
sheet [Birn et al. y 1997; Borovsky et ai, 1997]. 

The first row shows heating after steady state flux levels 
have been reached with constant Kp-\. Notice the strong day- 
side heating from the PEs with just a narrow band of heating 
around the nightside near L= 5. Heating from the PSEs is lim- 
ited to a strip around the dawnside along the outer boundary. 
This heating from the plasma sheet is so much weaker than the 
PE heating, it hardly changes the dawn sector of the summed 
plots from that of the PE results. It is interesting to note that 
most of the nightside is insignificantly heated by superther- 
mal electrons for this case. 

The second row shows results for steady state superthermal 
electron fluxes after constant Kp- 6. The PE heating is greatly 
diminished from the Kp - 1 values, with essentially no PE heat- 
ing on the nightside. By contrast, though, the PSE energy 
deposition has increased dramatically throughout the inner 
magnetosphere. Keep in mind that these values are propor- 
tional to the thermal plasma density, so the decrease at large 
radial distances is primarily due to decreases in this density. 
This is also true for the deposition in the afternoon sector from 
both populations. The dependence of P e on the thermal elec- 
tron temperature is negligible because the characteristic en- 
ergy of the superthermal electrons is much greater than T e . 

The bottom three rows of Plate 1 are results using the Kp 
history during the CRRES banded structure observations in late 
January 1991. The third row is during the final injection of 
plasma sheet particles when the peak ot the cloud is at local 
midnight, the fourth row is twelve hours later, and the final 
row is more than two days after this, chosen so the cloud peak 
is again at local midnight. In the third row, significant heal- 
ing is seen around the nightside due to PEs, in fact in all three 
times the nightside heating is more than the steady-state heat- 
ing when Kp is held at 1. This is mostly due to differences in 
the thermal plasma density between the two simulations. The 
PSE heating rates are less than the Kp=6 values, and are not so 
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close to the Earth, but are still substantial, particularly with a 
denser source population. Twelve hours later, the cloud is on 
the dayside and far less of an influence because the PE heating 
is so much stronger there. In the fifth row, the band is degrad- 
ing but still significantly contributing to heating the thermal 
plasma on the nightside, and is certainly still more important 
than PEs in this region. 

6.2, Flux Comparisons 

The heating rates in the previous section describe the com- 
posite influence of the superthermal electrons on the thermal 
plasma. It is also interesting to examine the combined distri- 
bution function in the superthermal energy range to understand 
how these bulk influences are generated. This will be done for 
the January 1991 results. As seen above, the Kp= 1 results 
have no PSEs in most of the simulation domain and the steady 
state Kp=6 case is somewhat unrealistic (shown as an extreme 
case). Therefore these results will concentrate on a compari- 
son of the distributions of PEs and PSEs during the CRRES 
banded structure observations. 

Figure 15 shows energy spectra for the two source popula- 
tions at various spatial locations for the third row of Plate I , 
during the final Kp spike and PSE capture. The spatial loca- 
tions are chosen to be in the nightside peak of the cloud (first 
row) and the nightside PE band (second rwo), and two analo- 
gous points on the dayside. Also, two equatorial pitch angles 
are shown: 90° and the loss cone boundary (LCB). Row 1 
clearly shows the PSE dominance of the distribution through- 
out the pitch angle range, with fluxes several orders of magni- 
tude larger than the PE fluxes at this location. In the narrow PE 
band on the nightside, however, PEs dominate the low-energy 


PA=90 PA=LCB 



Figure 15. Energy spectra for photoelectrons (solid line) 
and plasma sheet electrons (dotted line) at various pitch angles 
and spatial locations on January 27, 1991, at UT=1800. This 
time corresponds to the third row of Plate 1. 


range of the spectrum, even at 90°. On the dayside, however, 
the intersection shifts to higher energies because of the prox- 
imity of the PE source, reaching 400 eV near the LCB at L=4. 
In every erse, however, the PSEs dominate the high-energy 
part of the ;pectrum, often with positive gradients in the dis- 
tribution at the intersection and elsewhere. It can be seen that 
rotating the PSE results by 12 hours will make the energy of 
the PSE-PE intersection decrease in the dayside plots and will 
also make the PEs more comparable to the PSEs on the night- 
side. As time continues, the PSE fluxes will degrade (no more 
injections), allowing more of the PE distribution to dominate 
the combined flux function. 

The combined distribution function of superthermal elec- 
trons from these two sources is shown in Figure 16 throughout 
the CRRES observations at the 4 spatial locations discussed in 
Figure 15 (again at 90° and the LCB). Here the two popula- 
tions have been summed into a single distribution. In the first 
column, PShs form the nightside and PEs form the dayside dis- 
tributions, except for the high-energy bulge at 90° at L= 6. Af- 
ter this, the distribution becomes mixed, except in the first 
row, which ;s always dominated by PSEs. Notice the substan- 
tial number of spikes forming in the two L= 4 rows at high en- 
ergies. Th s is the banded structure forming inside of the 
AIfv£n boundary after the capture of the PSEs, as the magnetic 
drifts cause the electrons at these energies to have a drift period 
shorter than the corotation period, and they eventually lap the 
low-energy flectrons. This is not seen for most of the L= 6 re- 
sults because this is typically beyond the Alfven boundary. 
Only after extended quiet times would this radial distance de- 
velop the banded structure seen closer in. Also notice that the 
pitch angle distributions always show' a positive slope at the 
point of PE-PSE crossover. Furthermore, this slope is usually 
quite steep. 

7. Discussion 

By comb ning two time-dependent kinetic models of super- 
thermal electron transport, a global calculation of the super- 
thermal electron distribution function throughout the inner 
magnetosphere has been conducted. It has been shown that the 
energy range of validity for this combined model extended 
down to the superthermal-thermal intersection at a few eV 
(omitting transient field-aligned flows), allowing for the cal- 
culation of t ie entire distribution function and thus an accurate 
heating rate :o the thermal plasma. Because of the linearity of 
the formulae the source terms were separated to calculate the 
distributions from the various populations, namely pholoelec- 
trons (PEs) rnd plasma sheet electrons (PSEs). These distribu- 
tions were recombined to show the development of the total 
supertherma electron distribution function. 

It was sh iwn that convection, corotation, and Coulomb 
collisions an: the dominant processes in the formation of the 
PE distributi in function. Because the source is the dayside 
ionosphere, a source cone exists in the magnetosphere that 
fills in the cayside trapped zone. These electrons then propa- 
gate around he nightside, experiencing depletion due to colli- 
sions with the thermal plasma and flux tube relaxation from 
magnetospheric convection. As a result, only a small band in 
L shell is stil populated by PEs in the predawn sector, and this 
is mostly high-energy electrons in the trapped zone. This 
population creates a bulge in the PE distribution function on 
the morningsjde. Convection is also responsible for prevent- 
ing the PEs from building up completely on the dayside be- 
cause it sweeps away the particles from the trapped zone as 
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Figure 16. Time development of the combined-source superthermal electron distribution function at several 
spatial locations and pitch angles. 


they are scattering in from the source cone. This effect is rela- 
tively small during quiet times, but becomes a dominant influ- 
ence on the distribution function formation during disturbed 
conditions. It was shown that magnetic gradient-curvature 
drifts have essentially no effect on PEs because they are too 
low in energy to have a significant magnetic drift. 

PSEs are dominated by the interplay between the drift terms. 
Coulomb collisions only contribute to the formation ot the 
low-energy range of the distribution, and most of the high-en- 
ergy portion is not greatly affected by these interactions. Un- 
der steady geomagnetic conditions, a flow pattern forms that 
prevents the PSEs from existing for long periods in the inner 
magnetosphere. When the activity increases, PSEs push 
closer in toward the Earth. When the activity subsequently re- 
laxes, the PSEs closest to the Earth are then captured inside the 
Alfv6n boundary, where they will rotate and remain for several 
days. Magnetic drifts cause the high-energy electrons to su- 
percorotate, eventually lapping the lower-energy electrons and 
forming a banded structure in the energy spectra. 

The combination of these two distributions shows that PEs 
supply a low-energy population to the distribution function, 
but that PSEs dominate the high-energy part of the distribution 
function, particularly at large pitch angles. On the dayside, 
there is a change from a source cone distribution to a loss cone 
distribution as energy increases above the PE range. On the 
nightside, however, the fluxes at small pitch angles are rarely 
larger than those deep in the trapped zone. It is concluded that 
the addition of PEs does not greatly alter the nightside mean 
energy of the captured PSE cloud observed by the CRRES satel- 
lite. Because PEs exist in such a narrow radial band in the 
nightside, this cannot account for the less-than-adiabatic in- 
crease of T ± calculated from the lowest energy band observed 
by CRRES. Other mechanisms must be investigated to fully 
understand this phenomenon. 

A distinctive feature of most of the simulations presented in 
this study is an absence of superthermal electrons inside of 


about L - 3 on the nightside, prohibited by Coulomb losses 
(PEs) or corotation (PSEs). Only during prolonged periods of 
high geomagnetic activity can PSEs reach these L shells, and 
even then heating in this region is usually much less than on 
the dayside. An analysis of observed thermal plasma tempera- 
tures in this region for this effect is beyond the scope of this 
study, however. Conversely, these results show that PEs pump 
a substantial amount of energy into the thermal plasma on the 
dayside, particularly in the afternoon sector. These energy 
fluxes into the ionosphere typically approach 10 10 eVcm V 1 
in the plasmaspheric region, dropping to 10 8 and less outside 
the plasmapause. The PSEs, however, usually do not provide 
more than 10 8 eVcm'V 1 to the thermal electrons for the cho- 
sen boundary conditions. 

Note that while the combined model is capable of calculat- 
ing the precipitation effects of superthermal electrons into the 
ionosphere, this study did not include this process for the sake 
of brevity and clarity. However, it is expected that the inclu- 
sion of this population will substantially enhance the energy 
deposition rate from the PSEs because these secondary elec- 
trons will have a much lower mean energy and therefore more 
efficiently heat the thermal plasma. It is also expected that 
this population will increase the extent of nightside heating 
from the PEs. 

Another result from these combined spectra is that the dis- 
tribution often has positive slopes with energy or pitch angle, 
indicating that instabilities could arise. Because the present 
study has focused on global convection along and across field 
lines, an analysis of the stability of these distributions has 
not been conducted here. This will be the topic of a later study, 
including the feedback from internally-generated waves and the 
influences of externally-imposed waves. 
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